MD simulations were performed with the NAMD_2.12b1_ Linux-x86_64-multicore-CUDA suite of programs using ff14SB^[25] and bsc1^[26] as the protein and nucleic force field and TIP3P water molecules^[27] for water and Jong/Cheatham parameters^[28–31] for monovalent ion and Li/Merz^[32,33] (12-6 LJ nonbonded model) ion parameters for Mg²⁺ ion. The initial coordinates for the MD simulation were taken from the nuclear magnetic resonance (NMR) structure of CREB/CRE complex (PDB ID: 1DH3) as Fig. 1.^[19] Two different systems were employed, one with Mg²⁺ ion, and the other without Mg²⁺ ion. Both systems were prepared using the TLeap program in Amber16 package^[34] to resemble the cellular environment. The hexahydrated Mg²⁺ ion was bound in the cavity between the basic region of CREB and CRE. In the system, a solute complex was solvated in a period cuboid box, ensuring at least 12 angstrom water shell around the solute atom; the closeness parameter was 1.0 for solute atom and water. We started the simulation using a hexahydrated form of the ion, which may reduce the effect of force field artifact.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The crystal structure of CREB and CRE complex. The rose red ball refers to Mg2+ ion, and the red ball refers to water molecule.

Two systems were initially subject to energy minimization with harmonic restraints of on solute backbone atoms, and then, restraints would be relaxed during another energy minimization. After the energy minimization, two systems were heated to 310 K under canonical ensemble (NVT) during 620 ps simulation with harmonic restraints of on solute backbone atoms followed by 2 ns constant temperature (310 K) simulation. Below, five parallel MD simulations would be performed for each 250 ns with the temperature maintained at the set temperatures by Langevin dynamics and the Langevin piston Nose–Hoover method.^[35–37] The dynamic collision frequency was 1.0 ps, and the pressure was maintained at 1.0 atm using isotropic positional scaling with a pressure relaxation time of 2.0 ps.

To quantitatively examine the structural stability of CREB/CRE, we calculated the root-mean square deviation, the number and occupancy of hydrogen bonds, the width of the minor and major grooves, and the intra-base parameters of DNA.

RMSD^[38] can measure the deviation of a set of coordinates at time t_j to a reference set of coordinates, which are usually the coordinates of the given set of atom coordinates during a time of N_t. The RMSD is defined as 1 2 The RMSD values were calculated every 0.5 ps using the module cpptraj from Amber 16.

To understand thermodynamic stability, we calculated the number of hydrogen bonds every 0.5 ps using the visual molecular dynamics (VMD) plugin hydrogen bonds. The standard of hydrogen bond was set up as follows: the distance cutoff (acceptor to donor heavy atom) should be less than 3.5 Å; the angle cutoff (donor–hydrogen–acceptor angle) should be more than 140 degrees. To derive the average fluctuations in DNA conformations during the MD simulations, we have used the Curves+ program^[39] of Canal and Canion utilities.

The dynamic stability of CREB/CRE complex is elucidated by calculating the RMSD values of the backbone of protein and DNA with respect to their crystallographic coordinates. The entire simulation trajectory for 250 ns is calculated in MD simulations.^[38] The plot of these RMSDs as a function of time is provided in Fig. 2. The average value of RMSD for backbone atoms in the CREB/CRE with the Mg²⁺ ion is 2.9018 Å, whereas the average value of RMSD for backbone atoms in CREB/CRE without the Mg²⁺ ion is 3.3254 Å. The RMSD value plot shows that within 200–250 ns, the average RMSD value of the complex fluctuates slightly. It indicates that the complexes in the presence and absence of Mg²⁺ ion are stable, and Mg²⁺ ion has the same fluctuation during the simulation. At 250 ns, the system without the Mg²⁺ ion is 0.5 Å higher than the RMSD value of the system with Mg²⁺ ion, hence the system is more stable in the presence of the Mg²⁺ ion. The relatively high values and the irregular profile of the plots reflect the structural changes of highly flexible regions. Hence, it is reliable to perform the hydrogen bond occupies analysis from MD trajectory.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. RMSD analysis of trajectories of CREB/CRE complex as a function of simulation time. The case with the Mg2+ ion is indicated with black solid line. The case without the Mg2+ ion is denoted with red solid line.

It is important to recognize the main residues which play a key role in the binding between CREB and CRE. To understand the restricted nature of the CREB/CRE complex, it would be interesting to explore the hydrogen bond between DNA and adjacent protein residues. In the experiment, the authors analyzed the hydrogen bond formed in proteins. Hydrogen bonds are the most important specific interactions in biological recognition processes and they are one of the major driving forces for complex stability.^[24] Hydrogen bonds can be defined by using geometric or energy standards. In this study, we have used geometric criteria to define hydrogen bonds. We have analyzed the hydrogen bonds generated by CREB and CRE at the interface. Through all atom molecular dynamics, two important differences between the CREB/CRE complex in the presence and absence of the Mg²⁺ ion are identified.

First, it is found that the average numbers of hydrogen bonds for the CREB/CRE are 31.73 and 28.92 during the last 50 ns of simulation time, respectively. The decrease in the average number of hydrogen bonds for the CREB/CRE complex without the Mg²⁺ ion indicates a structural instability of the DNA–protein, which is consistent with our RMSD analysis. By counting the distribution of the number of hydrogen bonds in each trajectory, we can see an increase in the number of hydrogen bonds when Mg²⁺ ion is added to CREB/CRE complex from Fig. 3. We have calculated the contribution of the hydrogen bond to the interaction of CREB/CRE complex with an occupancy rate greater than 30% on the last 50 ns trajectory. The results are presented in Table 1. We determine that it is a significant difference in the hydrogen bond occupancy of CREB/CRE in the presence and absence of the Mg²⁺ ion. In the absence of Mg²⁺ ion, we can find that almost all hydrogen atoms which form hydrogen bonds with the base of DNA are from arginase (ARG). ARG is positively charged under physiological ph conditions, and it attracts hydrogen bonds with negatively charged nucleic acids. When the Mg²⁺ ion is added to the system, some of hydrogen bonds formed by the ARG are broken. Because the Mg²⁺ ion is strongly positive, and some neutral amino acids such as ASN51, ASN116, LYS117, LYS118, and SER113 form new hydrogen bonds. Moreover, the number of newly generated hydrogen bonds is significantly larger than the number of broken hydrogen bonds. For example, residue ARG99, conserved across the basic leucine zipper (bZIP) family, does not create hydrogen bonds with CRE when the Mg²⁺ ion is absent, but does participate in two hydrogen bonds (to C35 bases) when the Mg²⁺ ion is presented.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The distribution of the number of hydrogen bonds in each trajectory. The black line indicates the case with the Mg2+ ion, and the red line indicates the case without Mg2+ ion.

Table 1.

Table 1.

Table 1. Average occupancy of hydrogen bond between CREB and CRE in the presence and absence of the Mg2+ ion during MD trajectory. . Acceptor Donor Occupies Mg2+ Without Mg2+ DC32-OP2 ARG114-NE 0.82 0.72 DC32-OP2 ARG114-NH2 0.73 0.66 DT13-OP2 ARG52-NE 0.88 0.91 DT13-OP2 ARG52-NH2 0.75 0.75 DG12-OP2 ARG56-NE 0.44 0.43 DG12-N7 ARG59-NH1 0.51 0.47 DT29-OP2 SER58-OG 0.41 0.43 DG6-N7 ARG102-NH1 0.78 DG6-O6 ARG102-NH2 0.73 DG33-OP1 ARG111-NH2 0.60 DG33-OP2 ARG114-NH2 0.55 DG27-N7 ARG47-NH1 0.58 DG27-O6 ARG47-NH2 0.63 DG12-OP2 ARG56-NH2 0.56 DC35-OP2 ARG99-NH1 0.50 DC35-OP2 ARG99-NH2 0.45 ASN51-OD1 DC28-N4 0.59 ASN106-OD1 DC7-N4 0.89 DG12-O6 ARG59-NH2 0.78 DG9-OP2 LYS117-NZ 0.67 DA31-OP1 LYS118-NZ 0.61 DT8-OP2 SER113-OG 0.78 DG6-OP2 ARG102-NH1 0.44 DT34-OP2 ARG107-NE 0.42 DG27-O6 ARG47-NH2 0.40 DC11-OP2 ARG59-NE 0.81 DC11-OP2 ARG59-NH2 0.70

Table 1. Average occupancy of hydrogen bond between CREB and CRE in the presence and absence of the Mg2+ ion during MD trajectory. .

Second, through VMD, we can observe the specific location of the hydrogen bonds. An interesting phenomenon is illustrated. Not only the number of hydrogen bonds is increased from 30 to 40 upon the addition of the Mg²⁺ ion, but also in some cases their compositional atoms are changed. We mark the positions of generated hydrogen bond on the bases of the DNA, and the result is shown in Fig. 4. It could be found that the hydrogen bonds between DNA and protein are mostly located near the CRE-conserved -TGACGTCA- sequence. Probably because of the presence of Mg²⁺ ion, the distance between the CRE and the CREB is brought closer in spatial position. As shown in Fig. 5, the hydrogen bonds are produced between ARG59 and the phosphate group of C11, which are conserved across the bZIP protein family. However, the acceptor of the hydrogen bonds from ARG59 is changed from phosphate group of C11 to O6 of G12 once Mg²⁺ ion is presented at CREB/CRE complex. The positions of hydrogen binding are different, some are combined with the phosphate backbone of DNA, and others are bound to the base of DNA. The large number of newly formed hydrogen bonds and the changes of hydrogen binding positions, together, may contribute to the increase in the stability and binding affinity of CREB/CRE complex with Mg²⁺ ion.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. DNA sequence used for the MD simulations. The upper strand is numbered with 6–17 in the –3 sense and the lower strand is numbered with 27–38 in the –5 sense. Hydrogen bonds with Mg2+ ion and without the Mg2+ ion observed during the MD simulations are indicated by black dots and hollow dots, respectively.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The existence of the hydrogen bond between protein interface residue and DNA base in CREB/CRE complex. (a) The hydrogen bond of CREB/CRE complex without the Mg2+ ion (R59NH2-C11OP2, R59NE-C11OP2). (b) The hydrogen bond of CREB/CRE complex with the Mg2+ ion (R59NH2-G12O6). The rose red ball indicates the Mg2+ ion.

CREB is a cellular transcription factor. When activated, CREB protein recruits other transcriptional coactivators to bind to CRE promoter upstream region. The protein also has an Mg²⁺ ion that facilitates binding to DNA.^[40] As depicted in Fig. 6, at 250 ns, it produces a significant DNA deformation of both CREB/CRE complexes in the presence and absence of the Mg²⁺ ion. We can see that it is apparent that the CREB/CRE complex with the Mg²⁺ ion interacts more compactly than without the Mg²⁺ ion. To understand the specific binding nature of cation interaction with duplex DNA is very important. In this case, the hexahydrated Mg²⁺ ion bound waters can form hydrogen bonds with base edges and phosphate groups. We focus our attention on the width variation in the major and minor grooves, and the effects of cations on groove geometry.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Comparison of the structure of CREB/CRE complex in the presence and absence of the Mg2+ ion at 250 ns. The blue color indicates the complex with the Mg2+ ion, and the green color indicates the complex without the Mg2+ ion. The rose red ball denotes the Mg2+ ion.

In order to identify how the Mg²⁺ ion affects the binding of CREB to CRE, we have analyzed the average changes of the DNA major groove and minor groove width during the last 50 ns of simulation time. From Fig. 7(a), it is illustrated that the width of the minor grooves near the Mg²⁺ ion becomes large. From Fig. 7(b), it is shown that the width of the major grooves near the Mg²⁺ ion becomes small. The conserved TGACGTCA sequence provides an example of the cation-specific nature of major groove coordination. For this sequence, the Mg²⁺ ion is bound in the major groove near the center of the two base pairs G9P–G30P. We proposed a possible mechanism for the major groove narrowed when the Mg²⁺ ion is presented in Fig. 7(c). There are three amino acids near the major groove, LYS117, LYS62, and ARG59, which form a stable triangular structure. When the Mg²⁺ ion is added to CREB and CRE complex, the distances between the N atoms on the above three amino groups become closer to each other by electrostatic interaction than before. We can also observe that the distances between G9P-LYS117N and G30P-LYS62N get closer, and the bases G9 and G30 are located on the phosphate backbone of DNA. Therefore, it is believed that the changes in distance between them are the result of the electrostatic interactions, which further leads to a reduction in the width of the major grooves of DNA and an increase in the width of the minor grooves of DNA. As presented in Fig. 7(d), we can see that the width of the major groove is reduced from 17.1 Å to 15.6 Å in the presence of Mg²⁺ ion, which confirms our previous analysis.

Fig. 7.

Figure Option ViewDownloadNew Window

Fig. 7. The variation of the width of the major and minor groove in the DNA. (a) The difference in the major groove width of each base. (b) The difference in the minor groove width of each base; we take the CREB/CRE complex without the Mg2+ ion as control. (c) The change in the major groove width of base G9–G30. (d) The change of the distances between the N atoms of LYS117 and LYS62, LYS117 and ARG59, LYS62 and LYS117. The number in the brackets is the distance of CREB/CRE complex without the Mg2+ ion. The number above the brackets is the distance of CREB/CRE complex with the Mg2+ ion.

In order to compare the conformational changes of CRE, we have analyzed the average changes of intra-base pair parameters, buckle, shear, and opening during the last 50 ns. The calculated results are shown in Fig. 8. The base pair parameters have a certain change between the CREB/CRE complex in the presence and absence of the Mg²⁺ ion, and the rotation angle becomes larger. Because we only analyze the effect of an Mg²⁺ ion in the crystal structure on the CREB/CRE complex, the changes are slightly small. In view of the paucity of amino acid-base interactions, the recognition pattern of CREB/CRE complex and the Mg²⁺ ion depends on the sequence of conserved TGACGTCA. It is possible that the electrostatic interaction of Mg²⁺ ion promotes the binding of Mg²⁺ ion with CREB/CRE by appropriately pre-bending DNA. It also causes the minor grooves to widen and the major grooves to narrow at the interaction site.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. DNA intra-base pair parameters of the CREB/CRE complex. (a) The difference in open of each base. (b) The difference in buckle of each base. (c) The difference in shear of each base. The CREB/CRE complex without the Mg2+ ion is chosen as control.